function LE_val = LE(params, y0, T_total, dt, h)
    % 使用 RK4 计算最大 Lyapunov 指数 (LE) 
    % 输入：
    %   params  - 参数结构体（传递给模型函数）
    %   y0      - 初始状态向量
    %   T_total - 总模拟时间
    %   dt      - 重新归一化间隔
    %   h       - RK4 积分步长
    %
    % 输出：
    %   LE_val  - 最大 Lyapunov 指数的估计值
    
    % 初始扰动幅度
    delta0 = 1e-6;
    y1 = y0;
    % 在初始状态上加上一个微小扰动（此处在第一个分量上加扰动）
    y2 = y1 + [delta0; 0; 0];
    
    num_steps = floor(T_total / dt);
    sum_log = 0;
    time = 0;
    
    for i = 1:num_steps
        tspan = [time, time + dt];
        % 使用自写 RK4 积分器计算两个轨迹
        [~, sol1] = rk4(@(t,y) hybridNeuron(t, y, params), tspan, y1, h);
        [~, sol2] = rk4(@(t,y) hybridNeuron(t, y, params), tspan, y2, h);
        y1 = sol1(end,:)';
        y2 = sol2(end,:)';
    
        % 计算两个轨迹之间的欧氏距离
        d = norm(y2 - y1);
    
        % 累计局部指数 log(d/delta0)
        sum_log = sum_log + log(d / delta0);
    
        % 重新归一化扰动方向
        direction = (y2 - y1) / d;
        y2 = y1 + delta0 * direction;
    
        time = time + dt;
    end
    
    % 最大 Lyapunov 指数估计：总增长率除以总模拟时间
    LE_val = sum_log / T_total;
end
